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NOMENCLATURE 


A  area 

D  distance,  objective  function 
E  Eq.  (4) 

F  Eq.  (6) 

F  Eq.  (7) 
i  index 
J  index 

mean  distance  between  particles 

N  number  of  particles 

Pj^  permutation  of  the  index  array 

Arjj  distance  moved 

At  time  between  frames 
Vjj  matrix 

Vj^  mean  velocity  of  particles 
Aj  Lagrange  multiplier 
/Zj  Lagrange  multiplier 
p  weight  factor 
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1.  INTRODUCTION 


Particle  tracking,  as  in  particle  image  velocimetry  (PIV),  requires  the  identification 
of  each  particle  in  one  picture  with  another  in  a  second  picture.  The  second  picture 
shows  the  particles  after  a  short  amount  of  time  has  elapsed  so  that  the  particles  have 
moved  from  their  positions  in  the  first  picture.  The  identification  problem  is 
particularly  difficult  when  the  particles  must  be  regarded  as  indistinguishable,  i.e., 
individual  shapes  or  intensities  cannot  be  used  as  reliable  cues  of  identification.  A 
further  difficulty  arises  when  a  history  of  particle  trajectories,  i.e.  many  pictures  of  the 
])article  field,  is  not  available  as  further  tracking  information.  The  tracking  algorithm 
described  below  was  developed  for  specifically  these  circumstances.  A  summary  of  past 
attempts  at  solving  this  problem  is  given  in  Ref.  1  and  2. 


2.  THE  TRACKING  PROBLEM 


2.1  Minimum  Movement  Criterion 

Let  there  be  N  particles  present  in  each  of  pictures  1  and  2.  Let  |  Ar  |jj  represent 

the  inferred  distance  of  movement  due  to  identifying  particle  i  in  picture  1  with  particle 
j  ill  picture  2.  It  is  assumed  that  both  j)ictures  have  been  digitized  so  that  all  possible 
|)arti(le  movements  \  Ar  |jj  are  known  and  tabulated.  For  any  set  of  pairings  i,  j(i) 

chosen,  there  is  an  inferred  total  distance  of  movement 


I)  = 


(1) 


Since  the  particles  are  assumed  indistinguishable,  and  assumimg  a  relatively  short  time 
duration  between  the  two  pictures,  it  is  logical  to  choose  a  pairing  rule  j(i)  such  that 
the  total  inferred  movement  of  the  particles  is  a  minimum. 


I)  —  m  i  n 
•’k 


N 


X  ^  ^ij(i)  ■ 


(2) 


'I'his  problem,  by  its  nature,  is  discrete.  The  number  j  to  associate  with  each  i  must 
be  an  integer  lying  between  1  and  N.  'riiercfore,  the  minimization  to  be  accomplished  in 
(2)  docs  not,  at  least  outwardly,  lend  itself  to  solution  by  a  continuous  minimum 
stoking  algorithm  such  as  Newton-Raphson.  On  the  other  hand,  the  Newton-Raphson 
technique  is  very  i)owcrful  and  easily  programmed.  For  these  reasons,  it  would  be  nice 
to  convert  our  discrete  problem  to  an  equivalent  continuous  one. 
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2.2  Continuous  Formulation 


A  method  for  accomplishing  this  is  as  follows.  Introduce  a  matrix  Vjj  which, 

temporarily,  represents  for  each  particle  i  a  fractional  or  partial  association  with  each 
particle  j.  Thus,  each  element  Vjj  is  selected  from  the  continuum  of  values  between  0 

and  1,  such  that 

0<Vj.<l,  (3) 


Now  the  total  distance  D  may  be  replaced  by  [compare  with  (2)] 


E=  IIVylArlij. 

i  j 


(All  sums  go  from  1  to  N,  unless  otherwise  indicated.)  Of  course  in  the  final  analysis, 
tor  each  i  all  Vjj  should  be  zero  except  for  one,  corresponding  to  the  choice  j(i),  which 

will  have  value  unity.  However,  with  a  minimization  technique,  such  as 
Ncwton-Raphson,  with  continuously  varying  parameters,  it  would  be  permissible  to 
allow  values  V-  away  from  target  values  0  or  1,  provided  they  converge  in  some 

limiting  sense  to  0  or  1  eventually.  This  is  the  philosophy  behind  the  approach 
described  next. 

The  Newton-Raphson  technique  lends  itself  nicely  to  equality  constraints.  The  loose 
requirements  (3)  may  be  further  refined  to  a  useful  set  of  equalities. 


N 

£  V..  =  1.  i  =  1. N, 

j 

(5) 

N 


The  problem  is  to  minimize  E  subject  to  the  constraints  (5).  Introducing  Lagrange 
multipliers  Aj  and  /ij  we  obtain  a  modified  objective  function. 
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=  minimum  (6) 


M 

1 

i  j  i 

1  j  J 

j 

1 

for  unconstrained  minimization  with  respect  to  V-,  Aj  and  //j. 

However,  there  is  one  problem  with  this  approach.  There  is  nothing  in  it  that  forces 
values  Vjj  toward  either  0  or  1.  There  is,  it  turns  out,  a  rather  simple  way  to 

accomplish  this. 


2.3  An  Entropic  Approach 


An  image  restoration  technique  was  introduced  some  years  ago,  (Ref.  3)  that 
confines  its  outputs  to  lie  between  0  and  1.  It  does  this  by  adding  an  entropy  sum,  plus 
a  "complementary  entropy"  sum,  to  the  quantity  that  is  to  be  minimized.  We  can  do 
the  same  here.  The  new  objective  function  is 


N  N  N 

N 

N-1 

N 

I-ir* 

Ivir> 

1  j  1 

j 

j 

i 

N  N 

i  j 


The  last  double  sum  consists  of  entropy  VlnV  terms  and  "complementary  entropy"  (1  - 
V)  In  (1  -  V)  terms.  Weight  factor  p  is  at  the  user's  disposal  and  is  taken  up  below. 

We  now  seek  the  calculus-extremum  solution  to  (7),  and  see  if  it  can  be  forced 
toward  0  or  1,  as  claimed.  Setting  djOWy^  =  0  of  Eq.  (7)  yields  the  intriguing  result 


1  -f  exp  [  {\Ip)  (I  Ar|j.  +  X-  +  \ 


i,  j  =  1,...N. 


(8) 


The  action  of  weight  p  is  now  apparent.  If  we  make  p  very  small  then  the  exponent 
tends  toward  either  -t-oo  or  -w  forcing  Vjj  toward  either  0  or  1,  just  as  we  want. 
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There  are  four  double  sums,  or  terms,  in  G,  Eq.  (7).  The  size  of  p  determines  the 
relative  weight  of  the  last  term  for  the  minimization.  Because  we  want  all  the  emphasis 
to  be  upon  the  first  terms  our  objective  function  F,  Eq.  (6),  we  want  p  to  approach 
zero. 

The  particulars  of  the  Newton-Raphson  approach  that  solves  the  problem  (7),  (8) 
are  as  follows. 


3.  Newton-Raphson  Tracking  Algorithm 

The  aim  is  to  arrive  at  a  set  of  parameters  {A^},  {p.-^}  such  that  given  by 
representation  (8)  satisfies  the  equality  constraints  (5).  This  is  accomplished  as  follows. 

(a)  Fix  weight  p,  e.g.,  p  =  10.  p  must  be  positive. 

(b)  Start  with  a  trial  solution  {Aj},  {/ij}. 

(c)  Form  Vjj  by  Eq.  (8),  all  i,j. 

(d)  Form  sums  tj  =  ^Vjj,  Uj  =  bJ 

j  i 

(e)  If  all  tj  and  all  Uj  equal  1,  STOP. 

(f)  If  not,  set  Atj  =  1  -  tj,  Auj  =  1  -  Uj,  all  i,  j.  (9) 

(g)  Solve  for  changes  AAj,  A/ij,  all  i,  j,  in  2N  equations 

j 

i 

The  coefficients  of  {AAj},  {A/ij}  comprise  matrix  (11)  below 

(h)  Update  Aj  -»  Aj  +  AAj,  p-^-*  +  Ap-^ 
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(i)  Go  to  (c). 


Once  the  loop  STOPS  ,  at  step  (e),  weight  parameter  p  is  reduced  by  a  factor  1/2,  and 
step  (c)  starts  a  new  Newton-R^phson  problem,  using  a  previous  solution  set  {Aj},  {/ij} 

as  the  new  trial  solution.  Weight  factor  p  is  reduced,  problem  by  problem,  until  all  V- 

values  are  sufficiently  close  to  either  0  or  1.  This  solves  the  correspondence  problem  (2). 


3.1  Accuracy 

Can  we  be  assured  that  a  stationary  solution  to  (7)  is  a  minimum,  and  in 
particular,  the  global  minimum?  To  judge  this  effect,  take  the  second  partial  derivative 
of  form  (7).  This  yields  a  quantity 


/>ll/Vij+l/(l-V,j))  (10) 

Since  Vjj  must  be  positive,  by  representation  (8),  and  also  p  is  positive,  quantity  (10) 

must  be  positive.  Therefore  function  (7)  is  concave  downwards,  and  there  is  only  one 
stationary  solution,  a  minimum.  Hence,  the  output  of  algorithm  (9)  must  always  attain 
the  absolute  minimum  to  the  problem. 

This  was  confirmed  by  experimentation  with  the  algorithm.  For  small  dimensioned 
problems,  where  N  <  8,  the  absolute  minimum  solution  can  be  obtained  by  simply 
trying  out  all  possible  particle  pairings.  In  all  cases  tried,  where  particle  positions  were 
generated  with  uniform  randomness,  the  solution  obtained  by  algorithm  (9)  matched 
that  by  simple  permutation.  A  case  N  =  20  is  shown  in  Fig.  1. 


Figure  1.  A  case  of  N=20  particle  pairs.  Arrows  designate  connections  established  by 
algorithm  (9). 
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3.2  Efficiency 

The  objective  is  to  use  the  algorithm  on  very  large  problems,  where  N  is  the  order 
of  1000  or  more  particles.  On  the  other  hand,  step  (g)  of  algorithm  (9)  requires 
inversion  of  a  2N  x  2N  matrix,  and  this  is  done  repeatedly  during  enaction  of  the 
algorithm.  To  avoid  the  requirement  for  multiple  inversion  of  2000  x  2000  (and  more) 
matrices,  two  tactics  were  taken.  These  dramatically  speed  up  execution  of  the 
algorithm. 


(1)  The  Vjj  values  that  approach  1  do  so  at  different  rates.  Some  do  after  only  3  or 

4  iterations,  while  others  require  40  or  more.  Also,  once  a  is  close  to  1,  it  rarely 

backs  up  and  approaches  1  ultimately.  Hence,  to  avoid  the  need  for  carrying  along  this 
unnecessary  baggage,  the  algorithm  removes  from  the  problem  any  i,  j  pair  whose  V  is 
close  enough  to  1,  and  then  continues  the  problem  with  the  remaining  elements.  In  this 
way,  the  problem  rapidly  '•educes  in  dimensionality  after  but  a  few  iterations.  For 
larger  problems  of  N  >  50,  the  saving  in  CPU  time  due  to  this  tactic  is  greater  than 
9C 


(2)  The  inversion  subroutine,  that  is  used  in  step  (g)  contains  three  nested  DO 
loops.  This  makes  it  the  major  user  of  time  for  the  algorithm.  However,  the  matrix  that 
is  inverted  has  the  special  form 


0  W„  W,2  ...  W,N 

0  W2,  W22  ... 

0 


0  0  .  .  Wn,  Wn2 

i 

w„  W2,  ...  Wn,  0 

j 


w 


NN 


0 


W,2  W22 


■  ■  ■  "w  0 

j 


0 


(II) 


W,N  W2N  ...  Wnn  0  0 
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where  W..  e  V..(  1  -  V..).  Since  all  V-j  are  positive,  the  diagonal  elements  of  this 
J  J 

matrix  dominate. 

On  the  other  hand,  the  innermost  DO  loop  operates  on  off-diagonal  elements. 
Hence,  it  was  possible  to  replace  this  DO  by  a  single  operation,  effectively  changing  the 
three  nested  DO's  to  two.  It  should  be  emphasized  that  the  result  is  now  only  an 
approximate  inversion  of  the  matrix.  However,  since  the  Newton-Raphson  algorithm  is 
iterative  anyhow,  it  can  recover  from  the  inaccuracy.  In  practice,  it  does:  the  method 
reduces  execution  time  by  about  30  to  40%. 

Another  opportunity  for  time  saving  grows  out  of  the  off  diagonal  zeros  in  the 
upper  left  hand  quadrant  of  (11).  This  permits  the  inner  DO  loop  to  be  contracted  by 
about  a  factor  of  25%,  and  without  further  loss  of  accuracy  in  the  inversion. 


3.3  Empirical  Results 


We  generated  random  particle  configurations  for  the  two  frames,  and  then 
processed  them  by  use  of  the  tracking  algorithm  (9).  With  up  to  N  =  8  particles  in  each 
picture,  results  were  confirmed  by  permutationsralgorithm  that  searches  for  the 
absolute  minimum  in  D  by  exhaustive  search  of  all  solution  space.  For  values  of  N 
greater  than  8  quantitative  assessment  of  the  method  is  possible  based  on  computer 
experiments,  such  as  simulated  Brownian  motion.  There,  both  the  initial  and  final 
position  of  particles  are  known  and  the  accuracy  of  the  matching  by  the  algorithm  can 
be  quantified. 

Three  cases,  consisting  of  100,  250  and  1000  particles  each,  see  Figs.  2,  3  and  4, 
were  run.  The  particles  were  displaced  by  amounts  determined  by  a  random  number 
generator  and  the  algorithm  was  used  to  attempt  matching  of  the  particles.  The  results, 
in  CPU  seconds  on  a  Cray-2  running  UNICOS  5.0,  may  be  summarized  as  follows: 

TABLE  1. 


No  of 
particles 

Case  1 

Case  2 

Case  3 

100 

4.6 

4.1 

8.2 

250 

43.5 

27.6 

31.5 

1000 

1329.0* 

283.7 

237.3 

*  This  case  did  not  run  to  completion.  In  this  length  of  time  the  problem  had  been 
reduced  to  138  particles. 

**  Sample  cases  with  three  different  matrix  inversion  strategies  were  run.  Exact 
inversion  was  used  in  case  1.  In  case  2  inexact  inversion  was  used  until  the  number  of 
iterations  exceeded  9.  Case  3  employed  inexact  inversion.  In  each  case  three  samples 
consisting  of  100,  250  and  1000  particles  were  analyzed.  The  results  are  summarized  in 
Table  1  showing  the  CPU  seconds  utilized. 
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The  memory  requirements  for  the  test  cases  that  were  run  on  the  Cray-2  are; 


TABLE  2. 


100  particles: 
250  particles: 
1000  particles; 


158,345  words 

582.321  words 

8.102.321  words 


Each  word  consists  of  64  bits  (8  bytes). 


For  reasonable  expectation  for  a  solution  to  the  identification  problem,  several 
conditions  need  to  be  met: 

(1).  In  their  initial  distribution,  in  the  plane  of  observation,  the  particles  must  not 
be  packed  too  closely.  In  situations  where  particles  are  in  close  proximity  and  where 
there  is  not  a  well  defined  flow  direction,  such  as  for  example  in  Brownian  motion, 
crossing  of  trajectories  may  occur,  possibly  leading  to  erroneous  selection  of  terminal 
positions.  (2).  Condition  1  implies  that  the  maximum  allowable  displacement  between 
frames  should  be  limited.  It  has  been  found  that 


AtIVml  <0.5L,„  (12) 


is  a  reasonable  choice.  Here,  At  is  the  time  elapsed  between  frames  and  L^^  the  mean 
distance  between  particles  in  the  initial  distribution  and  v^  the  mean  velocity  of  the 
particles.  In  general,  the  mean  distance  between  particles  in  a  plane  can  be  expressed  as 


L  = 
m 


4A/7rN 


(13) 


where  is  the  mean  distance  and  A  is  the  area,  N  the  number  of  particles. 


(3).  Finally,  it  must  be  ascertained  that  there  is  no  temporal  variation  in  the 
number  of  particle  between  frames.  This  ensures  that  pairing  of  all  the  particles  is 
possible. 
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Figure  2.  A  case  with  N  =  100  particles.  Arrows  indicate  established  connections 
between  particles  in  the  two  frames. 
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Figure  3.  The  matchup  of  250  particle  pairs.  Arrows  indicate  matched  particles. 
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Figure  4.  The  case  with  N  =  1000  particles.  Arrows  indicate  matched  particle  pairs. 
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4.  CONCLUSIONS 


A  tracking  algorithm  has  been  developed  that  automatically  minimizes  total 
distance  of  movement  as  a  criterion  of  identification  of  particle  pairs.  The  algorithm 
always  produced  the  absolute  minimum  distance-pairings  as  its  solution.  Cases  of  up  to 
1000  particle  pairs  were  run  successfully,  and  with  small  CPU  time  requirement.  The 
algorithm  is  unsuitable  for  situations  where  the  volume  particle  density  is  high.  Effort 
is  under  way  to  extend  the  applicability  of  this  methodology  to  such  cases. 
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